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We propose and apply the finite-element discrete variable representation to express the nonequi- 
librium Green's function for strongly inhomogeneous quantum systems. This method is highly 
favorable against a general basis approach with regard to numerical complexity, memory re- 
sources, and computation time. Its flexibility also allows for an accurate representation of spa- 
tially extended hamiltonians, and thus opens the way towards a direct solution of the two-time 
Schwinger/Keldysh/Kadanoff-Baym equations on spatial grids, including e.g. the description of 
highly excited states in atoms. As first benchmarks, we compute and characterize, in Hartree- 
Fock and second Born approximation, the ground states of the He atom, the H2 molecule and the 
LiH molecule in one spatial dimension. Thereby, the ground-state/binding energies, densities and 
bond-lengths are compared with the direct solution of the time-dependent Schrodinger equation. 

PACS numbers: 02.70.Dh, 05.30.-d, 31.15.-p 



I. INTRODUCTION 

The two-time Schwinger/Keldysh/Kadanoff-Baym 
equations (SKKBE), e.g. [TJ |5J [3], allow for a quantum 
statistical analysis of nonequilibrium processes on 
microscopic footing. To great success, the one-particle 
nonequilibrium Green's function (NEGF) has been 
computed from the SKKBE for a variety of homoge- 
neous quantum systems, e.g. for nuclear matter |U |5J [B], 
the correlated electron gas [7], dense plasmas [3 |9j ITU] , 
or electron-hole plasmas [HI [12] [T3l HH [15], where 
different types of many-body approximations, by di- 
agram technique, have been included in a conserving 
manner. On the contrary, NEGFs, only in the recent 
decade, have challenged attention with respect to spatial 
inhomogeneity, exploring localized, finite and strongly 
correlated systems. Examples are electrons in atoms 
and small molecules \16\ [17], few-electron quantum 
dots [TBI EH] an d charge carriers in lattice and transport 
models such as strongly correlated Hubbard chains [20], 
molecular junctions [21 and quantum dot levels coupled 
to leads |22j. 

Although computational capabilities have been per- 
manently increasing in the recent past, NEGF calcu- 
lations remain a demanding task for finite systems: in 
particular — including electron-electron correlations in 
second Born approximation — highly excited states in 
atoms or time-dependent phenomena related to their ion- 
ization are generally difficult to access, and, only very few 
attempts have been made 231(51]. Also, the describability 
of specific correlation effects, such as two-electron reso- 
nances in dipole spectra [25]. remain so far unanswered in 
NEGF approaches as they require an accurate and exten- 
sive (large-scale) computation of the temporal evolution 
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following an intense external perturbation. 

All above-mentioned finite systems rely on general 
(semi-) analytic basis expansions of the nonequilibrium 
Green's function. Nevertheless, concerning the numeri- 
cal complexity associated with the NEGF, a basis repre- 
sentation reveals restricted capabilities. This affects, in 
particular, the spatial resolution, where a relatively small 
number of single-particle orbitals (typically ni, < 60 are 
feasible) are not appropriate to resolve the nonequilib- 
rium dynamics when, e.g. in atoms, occupations of highly 
excited states are non-negligible or ionization processes 
are involved. The same is the case when specialized basis 
sets constructed from Gauss-type or Slater-type orbitals 
or potential eigenstates are being used. To this end, ex- 
tremely large basis sets are needed and the system under 
investigation requires a large-scale treatment. 

Another option is provided by grid-based methods. 
However, for inhomogeneous systems, no direct solu- 
tion of the SKKBE with grid-based — and, in turn, 
finite-difference — methodologies has been performed so 
far, that systematically includes binary correlations and 
memory effects. This is due to the fact that numeri- 
cal grid methods allow for intuitive control but require 
small mesh spacings which become impractical for the 
compound structure of the two-space two-time Green's 
function: We note, that, in full three space dimensions 
(3D), the NEGF is an eight-dimensional complex func- 
tion. However, also in one spatial dimension, generally 
severe problems arise in the framework of spatially ex- 
tended hamiltonians, where particles may occupy broad 
domains in coordinate and/or momentum space. Thus, 
an alternative method for NEGF calculations is desir- 
able, to which one attributes more numerical flexibility 
and efficiency, and which has the ability to combine the 
advantages of non-existent grid and standard basis ap- 
proaches. 

In this paper, we develop such a computational method 
based on the finite-element discrete variable representa- 
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FIG. 1: In FE-DVR representation, the interval [0, Xo] is partitioned into n e finite elements [x 1 , x 1+1 ]. In each FE, n g generalized 
Gauss-Lobatto points (denoted x z m ) provide the basis for the construction of a local DVR basis set. rib denotes the dimensionality 
of the extended basis covering the whole interval. 



tion (FE-DVR), see Refs. [26j W\ and Sec. [HA] This 
method allows for an efficient solution of the two-time 
SKKBE for the one-particle Green's function, at least, in 
one spatial dimension. As general system, we, thereby, 
consider N interacting electrons, the non-relativistic 
hamiltonian of which reads 



h = t 
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with kinetic energy t, a possibly time-dependent poten- 
tial energy v, and the binary interactions described by ii. 
Except for their spin orientations, all electrons are con- 
sidered identical (in mass and charge), and, throughout 
the present work, we will use atomic units. 

The use of the FE-DVR provides analytical expressions 
for the kinetic and potential energy in NEGF calcula- 
tions for strongly inhomogeneous systems. But the main 
achievement of the present paper is the first realization 
of a grid-based NEGF approach together with a very effi- 
cient treatment of the binary interactions. Explicitl y, in- 
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stead of 0(nf) interaction matrix elements, see Sec 
our method requires only 0{n 2 ) elements, which, in 
dition, need not to be precomputed as before in a com- 
plicated manner. With regard to the SKKBE, the latter 
point directly leads to much simpler, semi-analytical for- 
mulas for the first- and second-order self-energies, which 
are independent of the explicit form of the interaction, 
Sec. |II C| With these remarkable scaring properties, the 
FE-DVR essentially reduces the numerical effort, such 
that considerably less storage memory and computing 
time is needed, and, hence, enables calculations on sig- 
nificantly larger, more extended systems than before. 

In Sec. |III| we demonstrate the power of the approach 
and compute the noncquilibrium Green's functions for 
the one-dimensional He atom and the neutral molecules 
H 2 and LiH (also in one spatial dimension) as function of 
the interatomic distance. In the course of this, we focus 
on the ground-state properties and compare the Hartree- 
Fock and the second Born approximation to the exact 
solution obtained from the full few-particle Schrodingcr 
equation. Ignoring the nuclear dynamics [i.e. in the Born- 
Oppenheimer scheme] , the exploration of noncquilibrium 
properties is straightforward within the formalism pre- 
sented. However, a detailed discussion is deferred to a 
forthcoming publication. 



II. FINITE-ELEMENT DISCRETE VARIABLE 
REPRESENTATION 



The finite-element discrete variable representation is a 
hybrid approach [28] which combines finite-element (FE) 
methods, i.e. spatial grids, and the discrete variable 
representation[29j (DVR). In a DVR basis, a similarity 
transformation allows us to replace matrix elements of 
local operators [of the coordinates] by their values on a 
relatively small numerical grid. The high degree of accu- 
racy of this procedure, widely used in quantum chemistry, 
manifests its usefulness in solving quantum mechanical 
problems [5U], 

For the direct solution of the few-particle time- 
dependent Schrodinger equation, e.g. Ref. |31j and 
references therein, the FE-DVR is highly effective — 
often in combination with time-dependent close coupling 
(TDCC) — due to the accuracy of the DVR on the one 
hand and the sparse character of FEs on the other. How- 
ever, these scaling properties, that enable a well paral- 
lelizable code [27], are less important for our application 
of the method. Instead, we focus on the benefits of the 
FE-DVR regarding the treatment of binary interactions 
and self-energies, which require the main computational 
expense within the framework of nonequilibrium Green's 
functions. 

The general idea how to combine FEs with the DVR to 
construct an extended basis is outlined in the following. 



Thereafter, in Sec. II B we discuss and give formulas for 
the relevant matrix elements of t, v and u, which are fi- 
nally used in the equations of motions for the one-particle 
Green's function, see Sec. Ill C 



Basis construction 



We divide the interval [0,x ], which is of physical and 
numerical relevance regarding hamiltonian ([lj and may 
be spatially extended, into n e finite elements with arbi- 
trary boundaries x° — < x 1 < x 2 < . . . < x n °~ 1 , x n " — 
xo, see Fig. [l] In each FE i, i.e. in [a; l ,a; l+1 ], we then 
construct a local DVR basis based on the generalized 
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FIG. 2: (Color online) Structure of a FE-DVR basis {XmO)} 
with n g = 4, i.e. 5 local DVR basis functions in each element. 
While the 'element' functions (solid) are defined in a single 
FE, the 'bridge' functions (dashed and dash-dotted lines) link 
two adjacent FEs. 



and are orthogonal with regard to the generalized Gauss- 
Lobatto quadrature, see Appendix. The 'bridge' func- 
tion (to = 0) in Eq. ([5]) extends over two adjacent el- 
ements [element i has overlap with element i + 1] and, 
hence, ensures communication between different grid do- 
mains i and i' and guarantees continuity of any expanded 
quantity or Green's function, cf. Sec. |II C| The 'clement' 
functions are zero at and outside the element bound- 
aries. Generally, in Eq. ^ and in the remainder of 
this paper, superscripts are labeling elements i ranging 
0, 1, 2, . . . , n e — 1, and subscripts are denoting local DVR 
indices to ranging 0, 1, 2, . . . , n s — 1, compare with Fig. [T] 
In the first (last) FE, the DVR basis function that is part 
of the left- (right-) hand bridge is removed, assuming the 
many-body wave function of system ([TJ to vanish out- 
side the interval [0, xq]. Hence, the total basis set has 
dimension 



Gauss-Lobatto points [ 



1 



and weights uu l m 



ci n = -{(x i+1 -x i )x m +(x i+1 + x i )} , 



Wm 

2 



(2) 



When using n g Legendre interpolating functions, the 
points x m {standard Gauss-Lobatto points) are defined 
as roots of the first derivative of Legendre polynomials 
P n (x) according to 



d.r 



and the associated weights are 



n g( n g + l)[ p n 3 0m)] 2 



(3) 



(4) 



In our approach, we use a DVR basis of equal size 
in each FE, see Fig. [2] The generalization to differ- 
ent numbers of basis functions per element is straight- 
forward, and only slightly alters the matrix elements in- 
volved, cf. Sec. |IIB[ The one-dimensional FE-DVR space 
is spanned by the set of orthonormal[32 functions 



/L-i(s) + /<T>) 



to = (bridge) 




else (element) 



with Lobatto shape functions [ 



r m (*) = n 



(5) 



(6) 



for x l < x < x l+1 and f^x) = for x < x % as well as 



x > x 



i+i 



which have the property fmi x m') = $u'&' 



n b = n e n g 



(7) 



We note, that, with our construction of the spatial grid, 
see Fig. [TJ the formula for the dimensionality slightly 
differs from Refs. [57]. Here, we do not, separately, 
define the size of the local DVR basis set, which would be 
rig+1, compare with Fig. [2] Moreover, a generalization to 
higher dimensions is possible by using a product ansatz 
for the coordinate functions \27\. 



B. Matrix elements of operators t, v, and u 

To perform NEGF calculations with respect to sys- 
tem (JT|), we need the matrix elements associated with 
the kinetic, potential and interaction energy operators 
referring to the chosen FE-DVR basis. Thereby, integra- 
tions over coordinate space are calculated by using the 
generalized Gauss-Lobatto (GGL) quadrature, and case 
differentiations arise from the fact that the basis func- 
tions Xm( x ) split into element and bridge functions. 

The potential-energy matrix — and the matrix of any 
other local operator — turns out to be diagonal with re- 
gard to elements i and local DVR basis indices to: 



<m 2 W = / dx X ^ 1 (x)v(x 1 t)xl? l2 (x) 
Jq 

— 3iii2 ~m 1 m 2 ^mi W ) 



(8) 



with 



v(x l m ,t) , to > 

4(*) = < v(xi a _ 1 ,t)wi_ l +v(xi +i ,t)wi +i 



m = 
(9) 



vi' w mm' 



Hence, Eq. (|8]) implies that the potential energy is simply 
represented by a vector of dimension n\>. 

The operator of the kinetic energy is non-local as it 
involves information of different points in physical space. 
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As a consequence, i^' 2 m2 is not diagonal. Particularly, 
any finite difference method applied to approximate the 



second derivative, cf. Eq. (10 1, must be carried out with 



great care, since the basis functions xln{ x ) given in FE- 



DVR representation are continuous but do not have con- 
tinuous derivatives at x % . Here, we follow the deriva- 
tion of Refs. [23] and [33] and obtain the block diagonal 
structure }27| of the kinetic-energy matrix as 



1 /" 2 



1 A - fil [711*1 „;* 
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r 



, mi > 0, m,2 > 

, mi = 0, 777,2 > 

, mi > 0, 772 2 = ' 

, 777l = r Tl2 = 



(10) 



where the quantity t l mim2 is connected to the first 
derivative [26 of the Lobatto shape functions via 



E 



dx 



dx 



(11) 



Eqs. ^ and ( 10 1 embody analytic formulas for the ki- 
netic and potential energy when a (finite-element) DVR 
basis is involved. 

The most attractive feature of the FE-DVR represen- 
tation is, that it can also be used together with the GGL 



quadrature to construct the matrix elements of the inter- 
action operator u which is non-local and of two-particle 
type. In general, the binary-interaction matrix elements 
(two-electron integrals) are carrying a set of four index- 
pairs (7,777) accounting for the two-particle character of 
the pair interaction, see first line of Eq. (12 1. However, 



in FE-DVR, using the separable form of the discretized 
interaction potential u{\x — x'\), see Eq. (14), we arrive at 



a very simple, semi-analytic expression for these matrix 
elements. This opens the way towards efficient NEGF 
calculations: 



x rxo 

• in 1 in j in .jii! ~ / dx I dx Xmi ( X ) Xm 3 i x ) u (\ x ~ x I) Xm 2 ( X ) Xm 4 ( x ) 

Jo Jo 

u ili 2 u i3i4 u rn±m 2 <Jm3m i u 'rn\m 2 ^ 



(12) 



with the remaining, full (kernel) matrix 

ra\vri 2 / , m 3 rmi m^i^m 2 m^ i \ ) 



being symmetric and of dimension 7-15 x rib- Here, the 
quantities a* m are the eigenvalues of the real matrix 

TJc \r, *\=u(\x l - x 1 ' , \) = a 13 /3?™? TO /3T"?, m ' 

(14) 



and /3^ m / are connected to the eigenvectors 0^ m > via 



l^rn'r 



, 777 > 



0" , = 



Pm'(n g -l) W n g -l + An'O 



W 



, 777 = 



+1 



(15) 



The key point is that the full rank representation (14 1 
enables us to factorize the two integrations in Eq. (12 1 



so that each integral can be separately performed by 
the use of the GGL quadrature. In turn, the evalua- 
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tion of the two-electron integrals Wm I 1 2 m 2 m 3m4 reduces to 
the computation of a simple matrix of dimension rib x rib, 
cf. Eq. pj). 



In summary, the effort of constructing the two-electron 
integral in FE-DVR representation becomes comparable 
to computing any single-electron matrix clement (such as 
the kinetic or potential energy) besides an additional but 
numerically elementary matrix diagonalization. More- 
over, Eq. ( 12 I is not only memory-friendly [the required 



memory scales with 0(n%) instead of 0(n^)\] but also 
permits a much more efficient evaluation of interac- 
tion contributions, especially, self-energy diagrams, see 
Sec. II C| This is due to the high degree of diagonality de- 
termined by the product of Kronecker deltas in Eq. ( 12 1. 



Also, it is favorable that the integrals do not depend on 
the explicit form of the pair interaction. 



C. Schwinger/Keldysh/Kadanoff-Baym equations 



The FE-DVR basis, as set up in Sec. II A| allows us to 
expand the one-particle noncquilibrium Green's function 
G(l,2) = -i {f c ip{i)^{l'))&\, with space-time argu- 
ments 1 = (%, t), 1' = (x',t') and spin omitted, as 

G & l ') = E E *™ (*) (*') <& 2 (t, t') . (16) 

The time-dependent coefficients g l m^ m2 {t,lf) are in gen- 
eral complex and vary on the complex Keldysh time- 
contour^ C. Further, G(l, 1') obeys the SKKBE[TJ |5J |3] 

{id t -H{l)}G(l,l')=S c (l-l') (17) 

d2S[G](l,2)G(2,l') , 



where H (S[G]) denotes the one-particle energy (self- 
energy), the time-integral is performed over C, and 
Eq. ( 17 1 is accompanied by its adjoint equation for the 
second time argument. Using Eq. (16 1, the SKKBEs 



transform into equations of motion for the matrix g [di- 
mension is rib x ni, with rib as defined in Eq. Q] and attain 
matrix form, where H, G and £ are to be replaced by 
their matrix components, 



G(l,2) 
H(l) 

S[G](1,2) 



a 1112 (t t') 



i,*i*a 
"77147772 



(*) 



_ J-«l*2 _|_ 7; *l l 2 

mi?Tl2 77117712 



(*) 



[g](t,t') 



(18) 
(19) 

(20) 



^77717712 \ U 1 L ) 



ycorr,2i22 

7771777,2 



(i,f) 



and all prod ucts are to be understood as matrix prod- 



ucts. In Eq. (19 1, h^ m2 (t) has the block-diagonal struc- 
ture impri nted by the kinetic energy, cf. Eq. (|10[). More- 



over, Eq. (|20| separates the self-energy £[g]^ 2 m2 (i, t') 



into Hartree-Fock (HF) and correlation parts, both of 
which are, generally, full [of dimension rib x rib] and func- 
tional of g. The Hartree-Fock self-energy S HF and the 
correlation self-energy S corr in second Born approxima- 
tion attain the form 



Ez-.iiis n i3i3 (f f+\ _ yiite „*a*i (f f+) 
u 'm 1 m 3 ym 3 m 3 \''T 1 J "mim 2 Sra 2 mi Ti 1 J 



S™ rr m T i2 (*> = E E 9Z 2 m 2 (*, <& 4 (*, - & 4 (*, <& 2 (*, *')} 9lZn 3 (f, t) < 3 m3 u 



(21) 



r ' l2M 4 , (22) 



237713 747774 



where a € {1,2} accounts for the spin-degeneracy, and t + 
indicates the limit t — > t + e>o from above on the contour 
C. Equilibrium initial correlations concerning S corr are 
treated in the mixed Green's function approach [36| [37} 
138] , where G and X have complex time-arguments i>o +i t 
with t S [— /3, 0] and /3 being the inverse temperature — for 
the full set of equations involved see e.g. Ref. [38 . 



The self-energy expressions ( |2l[ ) and ( |22| manifest very 
simple forms which arise from the subtle structure of the 
FE-DVR basis, compare with Refs. [T71[TH]. In the time- 



local HF part, Eq. (21 1, the Hartree term is completely 
diagonal [just as v in Eq. (|8l] requiring a single sum over 
the index pair (13,7713), and the exchange term involves 
only a product of two matrix elements. Note, that si- 
multaneous summations over i and m are equivalent to a 
single sum with rib elements! With this in mind, the eval- 
uation of the second Born self-energy, scales with 0(n\) 
implying only two summations per matrix element. This 
has to be compared with the general basis representation: 
there, two sums are required for each full vertex point 
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in the second-order diagrams and, additionally, a single 
sum is needed for the start- as well as for the end-point 
leading to an effort of 0(nf) in total for second-order 
self-energies. The simplification of this is a main result 
of the present paper and provides the basis for addressing 
new classes of problems, in particular, laser-atom inter- 
actions. 

In conclusion, using the FE-DVR representation in 
combination with the two-electron integrals u^*^ 3 ^ m4 
of Sec. |IIB| it is possible to rigorously reduce the com- 
putational complexity for inhomogeneous NEGF appli- 
cations, at least, in one spatial dimension. In par- 
ticular, with Eq. (22 1, the effort becomes compara- 
ble to that in lattice models, see e.g. [2U] [STJ |2"2"] . 
which, by construction, are computationally much sim- 
pler. Once the Green's function G(l,l') is computed 
from the matrix form of Eq. ( 17 1 , many observables 
such as the one-electron density n(x,t) = — iG(l,l + ), 
the time-dependent dipole moment (and in turn the 
polarizability[17.) or the total energy are accessible [39 . 



III. MODEL ATOMS AND MOLECULES 



Hartree-Fock 




n g [n b ] 


£g F [a.u.] 






4 [43] 
9 [98] 
14 [153] 


-2.22 

-2.224209- 
-2.2242096 


Second Born 


n g [n b ] 


# r-grid points 


P 2ndB ,, 1 
is gs [a. U.J 




14 [153] 
14 [153] 
14 [153] 
14 [153] 


101 
301 
601 
1001 


-2.23 

-2.2334 • ■ ■ 
-2.23341 ■ ■ 
-2.233418- 


TDSE (exact) 






IT.TDSE r „ 1 
ii gs [a. U.J 



-2.2382578 



TABLE I: Ground state energy E gB of the ID He atom 
(with fully converged decimal places) as computed from the 
Green's function in Hartree-Fock and second Born approxima- 
tion. The exact energy is obtained from the time-dependent 
Schrodinger equation (TDSE). 153 FE-DVR basis functions 
[at n e = 11] are adequate to reach the HF-limit and, thus, 
convergence with respect to the basis size. In second Born ap- 
proximation, about 600 points in imaginary time are needed 
for convergence in the fifth decimal place. 



In this section, we apply the FE-DVR representation, 



Eq. ( 16 1, to compute the nonequilibrium Green's function 
for atomic and molecular few-electron model systems. 
As atomic example, we discuss the one-dimensional he- 
lium atom (ID He), e.g. 001 113 EH ESJ HE] , which rep- 
resents the most elementary closed-shell system. This 
model of the 3D helium atom has been studied since the 
1970s and is known to reliably provide the qualitative 
features of the single- and double-ionization dynamics in 
intense laser fields [45] including the knee structure |44j . 
Moreover, it is still actively considered, e.g. EE] [47], as 
it serves as a fundamental 'testing ground' for multi- 
electron calculations. This issue is due to the presence of 
strong electron-electron (e-e) correlations which require 
a treatment beyond mean-field (HF) theories. In addi- 
tion to He, we discuss two molecular models with two and 
four electrons, respectively: The hydrogen molecule (H2), 
e.g. |S] EES IMS ED EH, and lithium hydride (LiH)[53]— 
again in one spatial dimension. The reason why we fo- 
cus on these atomic and molecular systems is twofold: 
(i) the long-range character of the ionic Coulomb po- 
tential (enhanced in ID!) proves the vital necessity for 
extended basis sets for the construction of which the 
FE-DVR is indeed well suitable, and (ii) the possible 
comparison to exact solutions, obtained from the time- 
dependent Schrodinger equation (TDSE), allows us to 
verify the quality of the involved many-body approxima- 
tions. Also, in the present paper, we restrict the NEGF 
calculations to the ground states. 

In hamiltonian ([!]), the helium atom is modeled by 
using v(x) — —Z[{x — so/2) 2 + l] -1 / 2 as regularized po- 
tential, where the atomic number is Z = 2. Thereby, 
the xo/2-shift ensures that the nucleus is situated in the 
center of the discretized interval [0, xq\. For the hydro- 



gen molecule and lithium hydride the coordinate is taken 
along the bond axis such that the potential is given by 
v d (x) = -ZiUx - (x + d)/2) 2 + I]" 1 / 2 _ z 2 [(x - (x - 
d) /2) 2 + l] -1 / 2 , where d denotes the interatomic distance, 
Zi=l and Z2 = 1 for the hydrogen and Z2 = 3 for the 
lithium atom. Principally, the regularization parameters 
[here, 1 for H and Li] can be adjusted to match the differ- 
ence of the ionization potentials of the individual model 
atoms to the 3D atoms, see Ref. [53]. Furthermore, for 
all three systems, a soft-core Coulombic e-e pair potential 

has been applied: u(\x — x'\) = Ux — x') 2 + l] 

For the ID helium atom, we have used 11 finite- 
elements within an interval of xq — 50 a.u. length. Some 
smaller FEs have, thereby, been arranged around xo/2 to 
ascertain larger numerical precision in the central region. 
Further, the number of local DVR basis functions n„ + 1 
has been varied between 5 and 20 to obtain convergence 
of the ground-state energy E gs , and, in Eqs. (21 1 and 
(22 1, the spin-degeneracy factor was set to a = 2 leading 



to the singlet state. 

For the Hartree-Fock approximation, the conver- 
gence of the He ground-state energy — at the fixed FE- 
configuration — is shown in Table [i] with regard to the 
basis size. At n g = 14 , corresponding to 153 basis func- 
tions in total, we obtain the HF-limit with more than 
six decimal places precision and, consequently, sufficient 
convergence with respect to the basis dimension. For 
the second Born approximation, we used the same FE- 
DVR set-up. However, due to the grand-canonical aver- 
aging involved in G(l, 1'), see definition in Sec. II C| the 
ground-state [equilibrium] Green's function has an addi- 
tional imaginary time argument r = t — t' € [— i/3, 0] C C. 
This has been discretized using a uniform power mesh, 
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FIG. 3: (Color online) One-electron ground-state density n(x) 
of the one-dimensional He atom, the H2 and the LiH molecule 
in Hartree-Fock (dashed) and second Born approximation 
(solid). The exact density (dotted) is obtained from the time- 
dependent Schrodinger equation (TDSE) with imaginary time 
propagation. The colored dots show the equilibrium positions 
of the ions separated by the bond-length db, Table [IT] and the 
gray curves indicate the associated potentials v(x) [in a.u. (left 
ordinate) but scaled by factor 0.35 and shifted]. 



for details see e.g. Refs. [HOGIS], and to ensure the zero- 
temperature limit, i.e. the ground state, we set [3 = 100. 
We note that, in the HF case, this grid is redundant as 
E HF (i, t') is local in time, cf. Eq. (21 1. For the second 



Born calculation, this implies, though, checking conver- 
gence with respect to a second parameter: the number 
of r-grid points, see Table [i] In 2 nd Born approxima- 
tion, the helium ground-state energy converges towards 
-2.2334 a.u., which is 0.0092 a.u. lower than the HF ref- 
erence value, and a comparable accuracy is obtained by 
using more than 600 time-grid points. With a deviation 
of less than 0.005 a.u., it comes close to the exact ground 
stated (-2.2383 a.u.), which follows from the TDSE. 
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4: (Color online) Hydrogen-hydrogen binding energy 
b as function of the interatomic distance d for the case 
of the binding singlet (| ][)) and anti-binding triplet, spin- 
polarized (|TT)) system. While the triplet system is less af- 
fected by correlations (see inset figure), the binding-energy 
curve of the singlet state is essentially improved against HF 
in the 2 nd Born approximation. The triangles mark TDSE 
results from Ref. |48| . The exact dissociation threshold is 
indicated by the horizontal line. For the values of the bond- 
lengths see Table [H] 



The one-electron ground-state density for the ID He 
atom is obtained from n{x) = — 1 G(x, t; x, t')| T ^n+io- 
and is displayed in the bottom graph of Fig. [3] The dif- 
ferences of both approximate results (dashed/solid line) 
and the exact density (dotted line) are most dominant 
within a small range of 0.5 a.u. around the ion position. 
As for the total energies, the second Born density im- 
proves the HF result and is relatively close to the exact 
density profile. 

For the hydrogenic and the lithium hydride system, 
the electron ground-state energy changes with distance d 
between the atomic nuclei. Hence, whether [or whether 
not] the individual atoms combine into molecules, de- 
pends on the H-H (Li-H) binding energy E^(d) which 
is the electron ground-state energy plus the interatomic 
repulsion[55] (Z\Z-i)ld. The computed binding-energy 
curves for H2 are displayed in Fig. [4] where a FE-DVR 
set-up similar to the helium case has led to convergent 
results. The singlet state |TI), again with a = 2 in the 



Bond-length db 




HF 


2 nd Born 


exact 




H 2 


1.9925 


2.0561 


2.151 




LiH 


3.3860 


3.5053 


3.6 •■ 


Binding energy Eh 




HF 


2 nd Born 


exact 




H 2 
LiH 


-1.3531 
-4.8534 


-1.3740 
-4.8886 


-1.391 
-4.91- 



TABLE II: Computed equilibrium bond-lengths db and cor- 
responding binding energies E gB (dh) + Z\Z%jd\, of the one- 
dimensional H2 and LiH model [all quantities in a.u.]. While 
the Hartree-Fock and second Born values are Green's function 
results, the exact values are obtained from the full solution of 
the few-particle TDSE. 



self-energies, is binding showing a minimum at a distance 
rfb in the exact result (dotted curve) and a well denned 
dissociation threshold (horizontal line). Also, the HF 
(dashed line) and the second Born approximation (solid 
line) confirm a substantial hydrogen-hydrogen binding, 
where 2 nd -Born correlations lead to a larger binding en- 
ergy that indicates an essential improvement of about 
60% in the HF energy discrepancy. However, for the 
singlet state, neither the Hartree-Fock nor the second 
Born approximation can accurately resolve the dissocia- 
tion threshold at —1.3396 a.u.. This is due to the fact 
that the closed-shell H2 molecule dissociates into open- 
shell fragments — single hydrogen atoms. Such a transi- 
tion can not be described within the semi-local (spin- 
restricted) approximations involved in the NEGF. We 
note that, the same problem is encountered in time- 
dependent density functional theory (TDDFT) using ex- 
act exchange |56j. Nevertheless, the different equilibrium 
positions of the nuclei (bond-lengths dh,) and the ground- 
state energies are not affected by this failure of the many- 
body approximations, see Table |ll| Overall, it turns out 
that correlations cause larger bond-lengths due to the 
lower electronic ground-state energy. 

We, in addition, have performed calculations for the 
spin-polarized or triplet H 2 system, \\\) with a = 1. The 
respective binding-energy curves in Fig. [4] show that, in 
contrast to the singlet state, as expected, it does not 
undergo molecular binding but behaves correctly in the 
limit d — ► 00. In particular, for all interatomic distances, 
the exact binding energy is well approximated within HF. 
Correlations generally improve the results (see inset of 
Fig. [4]) but play a minor role. This is typical for spin- 
polarized systems. 

The approximate and exact one-electron ground-state 
density for the H 2 singlet is shown in Fig. [3] (center 
graph) with respect to the corresponding equilibrium 
bond-lengths. Thereby, the gray curves illustrate the 
ionic potentials Vd b (x) on the bond axis. The exact 
density profile indicates onset of electron localization on 
the individual hydrogen atoms. This is not captured in 
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123456789 
interatomic distance d [a.u.] 



FIG. 5: (Color online) a) Li-H binding energy E^ lH as func- 
tion of the interatomic distance d. For the specific bond- 
lengths see Table [H] The compound dissociates at a threshold 
of . . . a.u.. For comparison, the dash-dotted line shows the 
binding-energy curve for the three-dimensional molecule [57] • 
b) Most relevant natural orbitals <f>i(x) for the LiH ground 
state at equilibrium bond-length as obtained from the TDSE, 
the HF and second Born Green's function [weighted by their 
occupation m]. 



HF and 2 Born approximation, which both lead to a 
smooth profile, but the trend towards a lower density 
between the nuclei becomes obvious. In particular, the 
ionic potential with the second-Born value of dh [com- 
pare also the dots in Fig. [3] is in good agreement with 
the TDSE result. 

The four-electron molecule, lithium hydride, serves 
as a simple example for the hetero-atomic dissociation. 
However, LiH, just like molecular hydrogen, dissociates 
into open-shell components — Li(3e) and H(le). Thus, in 
Fig. a), we obtain a similar behavior of E^ lH in HF 
and 2 nd Born approximation compared to H 2 against in- 
teratomic distance. Within the calculations, we used a 
basis consisting of 15 non-equidistant elements and up to 
15 local DVR basis functions for large values of d. For 
LiH, the inclusion of e-e correlations, improves the results 
such that the minimum in the binding energy becomes 
situated below the exact dissociation threshold. This is 
not realized in HF approximation. For reference, we, in 
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Fig. [5] a), also included the binding-energy curve for the 
three-dimensional counterpart [57] (dash-dotted line) of 
the ID model, which possesses a stronger Li-H bond at 
comparable internulcar distance. However, we note, that 
bond-lengths and binding energies are very sensitive to 
the softening parameters used in the ion and Coulomb 
potentials, e.g. Ref. [55]. For the specific values of and 
di, for lithium hydride, see Table [IT) The one-electron den- 
sity of LiH is plotted in the top graph of Fig. [3] where the 
lithium (hydrogen) atom is situated at negative (positive) 
x-positions, cf. the ion potentials Vd(x). In all considered 
cases, the density shows a clear minimum between the 
nuclei, and correlations mainly alter the density around 
the hydrogen atom. In particular, we highlight, that the 
second Born ground-state density is in surprisingly good 
agreement with the exact result. 

In order to identify the intra-molecular electronic 
structure more closely, we, in addition, have computed 
the most relevant natural orbitals (NO) for the ID LiH 
molecule, see Fig. [5] b). The natural orbitals 4>i(x), 
i = 0, 1, . . . , rib — 1, are obtained from the eigenvalue 
problem 



dx' p(x,x')(j) i {x') = Uiipi(x) 



(23) 



with density matrix p(x, x') = — i G(l, l')r^o+io- an d oc- 
cupations rii £ [0, 1]. For HF ground states, with f3 — > 00, 
we have rii = 1 for i = 0, 1, . . . , N — 1 and zero otherwise, 
where N is the number of electrons with the same spin 
projection. Hence, there are two fully occupied orbitals 
for the case of lithium hydride in Hartree-Fock approxi- 
mation, see the NO <fio and in Fig. [5]b). Correlation 
effects generally lead to occupations of more than two 
orbitals, cf. §2 in second Born approximation, compare 
with the exact result, and note that the orbitals have 
been scaled by y/nl. On the contrary, the two core elec- 
trons at the lithium atom, according to the localized NO 
4>q, are very little affected by correlations, which is re- 
vealed by the HF-, 2 nd -Born- and the TDSE-curves lying 
almost on top of each other. Further, the second natural 
orbital <p\ — with about 95-98% occupation and a node 
near the lithium atom — is shared between the nuclei and 
extends over several bond-lengths. In correspondence to 
the one-electron density in Fig. [3j the exact <f>i is well 
approached by the second Born approximation, which 
shows the correct trend in the central bond region. Also, 
the third NOs (f>2 are similar in shape. However, the de- 
viation in their occupations, is mainly responsible for the 
differences of the 2 nd Born to the exact result. Finally, 
all other natural orbitals (those which are not shown) are 
occupied by less than 1%. 



IV. CONCLUSION 

In this work, we have applied the finite-clement dis- 
crete variable representation (FE-DVR) to expand the 



one-particle nonequilibrium Green's function with re- 
spect to the two [one-dimensional] spatial coordinates. 
This procedure is highly favorable against a general basis 
representation: (i) conceptionally, it allows for an opti- 
mal and flexible combination of grid and basis methods, 
(ii), with respect to the NEGF of finite systems, a di- 
rect solution of the SKKBE within a grid-based hybrid 
approach becomes possible by (iii) an essentially simpli- 
fied treatment of all binary interactions. The latter point 
includes the description of particle-particle correlations, 
where the second-order Born self-energy in Sec. |II C| as 
the most basic model of correlations, attains a compara- 
bly simple form induced by the high degree of diagonality 
involved in the two-electron integrals, Eq. ( 12 1 , expressed 
in the FE-DVR picture. Also, due to the discretization 
in coordinate space, it is straightforward to change the 
one-particle potential v(x) or the specific form of the pair 
interaction u(\x — x'\). This is in striking contrast to a 
general basis, where to some extent enormous, extra nu- 
merical effort is required if the matrix elements and/or 
two-electron integrals are not analytically accessible and 
have to be precomputed. This, completely, drops out in 
the present approach. 

In summary, the developed method enables better per- 
formance with relation to larger accuracy and spatial res- 
olution, but at crucially lower numerical cost — regarding 
storage memory and computing time. In particular, this 
also holds true when spatially extended hamiltonians are 
being considered, as shown in Sec. |III| In turn, applying 
the FE-DVR, larger basis dimensions with a guide num- 
ber of Ub « 500 — 1000 become feasible, which implies 
an enhancement of more than one order of magnitude 
compared to a general basis approach. 

For illustration purposes, we have computed the 
nonequilibrium Green's function for simple but bench- 
marking atomic and molecular models: The helium atom 
and the linear molecules H2 and LiH in one spatial dimen- 
sion. Especially for the molecular systems, where two 
(four) electrons are shared between the nuclei, the en- 
hanced electron collision rate in one dimension makes it 
attractive to investigate electron-electron correlation ef- 
fects in second Born approximation. Indeed, with respect 
to inhomogeneous and finite systems, only few compar- 
isons of NEGF findings to exact many-body results are 
available. In our comparisons, we restricted ourselves 
to two- and four-electron models as the full solution of 
the TDSE becomes impractical for more than four elec- 
trons. In the present examples, it turns out, that the 2 nd 
Born approximation is well capable to catch the main 
ground-state features of the considered models. Thus, 
the presented analysis affirmatively contributes to the 
assessment of the applicability of NEGFs to atomic and 
molecular systems. 

Of course, the FE-DVR approach enables calculations 
also with larger particle numbers. Depending on the sys- 
tem, multi-electron ensembles [in one spatial dimension] 
with up to N < 20 turn out to be feasible [58 . Partic- 
ularly, we note that, with this grid-based method ad- 
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equate spatial resolution over a range of several hun- 
dred atomic units becomes, for the first time, available 
in NEGF approaches to strongly inhomogeneous quan- 
tum systems. The good performance is thereby not lim- 
ited to the second Born approximation. The method 
also allows for more complicated self-energies including 
GW or T-matrix on spatial grids. Moreover, the at- 
tractive scaling behaviors of the FE-DVR fully survive 
in nonequilibrium situations and, thus, provide essential 
impact for the efficient solution of the two-time SKKBE 
for atomic and molecular systems. Explicit results of the 
time-evolution in second Born approximation, including 
transitions to few-electron resonance states [25] located 
energetically above the one-electron excitations, will be 
the subject of a forthcoming publication. 
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APPENDIX: GENERALIZED GAUSS-LOBATTO 
INTEGRATION 

In numerical analysis, the generalized Gauss-Lobatto 
(GGL) scheme is a special quadrature rule which approx- 



imates a definite integral of a function g{x) as 

/ dxg(x) = ^2 dxg(x) 9{x % m ) w m > 

Jo i Jx% i m=0 

(A.l) 

where the specified points x l m and weights w l m are associ- 
ated with sub-domains [x l , x l+1 ] or finite elements i of the 
integration, see definition in Sec. |II A| For an arbitrary 
segmentation of the total domain [0,xo], the approxima- 
tion becomes exact in the limit n g — > oo. Moreover, from 
the GGL integration it follows, that the Lobatto shape 
functions are orthogonal in the sense of the quadrature 
rule: 

/ dx f rn (x) f m ,(x) 6jj> ^ ^ fm( X m) fm'( X m) W fh 
J fa 

= hi 1 S nirn> W % m ■ (A. 2) 
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